cd5 <- read.csv("HagPar-JLC-data.csv")

#############################################
##TABLE 2 (Results for Hypotheses 1 and 1A)##
#############################################

##Column 1
reg1a1 <- glm(fcomp~0+as.factor(postjud)+log(amount_use)+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg1a1)

##Column 2
reg1 <- glm(fcomp~0+as.factor(postjud)+log(amount_use)+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+govteff+rule_law+impunity+
              cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg1)

##Column 3
creg1a2 <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*revenue_world+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg1a2)

##Column 4
creg1a <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*revenue_world+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
                govteff+rule_law+impunity+cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg1a)


##################################
##FIGURE 1: Brazil and Guatemala##
##################################


#Data to use for plot
plotdat <- cd5[cd5$pnp==1,]

#Dummy variables for country so that we can take the "average" of this factor variable
plotdat$hon <- ifelse(plotdat$country=="Honduras",1,0)
plotdat$sur <- ifelse(plotdat$country=="Suriname",1,0)
plotdat$per <- ifelse(plotdat$country=="Peru",1,0)
plotdat$ven <- ifelse(plotdat$country=="Venezuela",1,0)
plotdat$col <- ifelse(plotdat$country=="Colombia",1,0)
plotdat$nic <- ifelse(plotdat$country=="Nicaragua",1,0)
plotdat$gua <- ifelse(plotdat$country=="Guatemala",1,0)
plotdat$els <- ifelse(plotdat$country=="El Salvador",1,0)
plotdat$chl <- ifelse(plotdat$country=="Chile",1,0)
plotdat$par <- ifelse(plotdat$country=="Paraguay",1,0)
plotdat$ecu <- ifelse(plotdat$country=="Ecuador",1,0)
plotdat$bol <- ifelse(plotdat$country=="Bolivia",1,0)
plotdat$pan <- ifelse(plotdat$country=="Panama",1,0)
plotdat$bra <- ifelse(plotdat$country=="Brazil",1,0)
plotdat$mex <- ifelse(plotdat$country=="Mexico",1,0)
plotdat$uru <- ifelse(plotdat$country=="Uruguay",1,0)
plotdat$arg <- ifelse(plotdat$country=="Argentina",1,0)
plotdat$cra <- ifelse(plotdat$country=="Costa Rica",1,0)
plotdat$dmr <- ifelse(plotdat$country=="Dominican Republic",1,0)
plotdat$hti <- ifelse(plotdat$country=="Haiti",1,0)
plotdat$tat <- ifelse(plotdat$country=="Trinidad and Tobago",1,0)

#logged award size
plotdat$lamuse <- log(plotdat$amount_use)

#drop incomplete observations
plotdat$keep <- ifelse(is.na(plotdat$amount_use==F),0,ifelse(is.na(plotdat$govteff==F),0,1))
plotdat2 <- plotdat[plotdat$keep==1,]

#model for plot (same as column 4 in Table 2)
plotmodel <- glm(fcomp~0+as.factor(postjud)+lamuse*revenue_world+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
                   govteff+rule_law+impunity+cumj_cy+hon+per+ven+col+sur+arg+cra+dmr+hti+tat+nic+
                   gua+els+chl+par+ecu+bol+pan+bra+mex,family="binomial"(link="cloglog"),data=plotdat2)

#plotting award sizes for predicted probability 
lmu <- seq(6.5,17,1.5)

#prediction model for Brazil (panel a)
p_bra <- predict(plotmodel,
                 data.frame(lamuse = lmu, 
                            revenue_world=rep(32.36,8),
                            postjud = factor(rep("5",8),levels=levels(factor(plotdat2$postjud))),
                            hon=rep(0,8),
                            per=rep(0,8),
                            ven=rep(0,8),
                            col=rep(0,8),
                            gua=rep(0,8),
                            els=rep(0,8),
                            chl=rep(0,8),
                            par=rep(0,8),
                            ecu=rep(0,8),
                            bol=rep(0,8),
                            pan=rep(0,8),
                            bra=rep(1,8),
                            mex=rep(0,8),
                            arg=rep(0,8),
                            sur=rep(0,8),
                            cra=rep(0,8),
                            dmr=rep(0,8),
                            nic=rep(0,8),
                            tat=rep(0,8),
                            hti=rep(0,8),
                            pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="Brazil"]),8), 
                            accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="Brazil"]),8),
                            npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="Brazil"]),8),
                            v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="Brazil"]),8), 
                            v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="Brazil"]),8),
                            rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="Brazil"]),8),
                            impunity=rep(mean(plotdat2$impunity[plotdat2$country=="Brazil"]),8), 
                            cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="Brazil"]),8),
                            govteff=rep(mean(plotdat2$govteff[plotdat2$country=="Brazil"]),8)),
                 type="response",na.action=na.pass)

#prediction model for Guatemala (panel b)
p_gua <- predict(plotmodel,
                 data.frame(lamuse = lmu, 
                            revenue_world=rep(11.51,8),
                            postjud = factor(rep("5",8),levels=levels(factor(plotdat2$postjud))),
                            hon=rep(0,8),
                            per=rep(0,8),
                            ven=rep(0,8),
                            col=rep(0,8),
                            gua=rep(1,8),
                            els=rep(0,8),
                            chl=rep(0,8),
                            par=rep(0,8),
                            ecu=rep(0,8),
                            bol=rep(0,8),
                            pan=rep(0,8),
                            bra=rep(0,8),
                            mex=rep(0,8),
                            arg=rep(0,8),
                            sur=rep(0,8),
                            cra=rep(0,8),
                            dmr=rep(0,8),
                            nic=rep(0,8),
                            tat=rep(0,8),
                            hti=rep(0,8),
                            pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="Guatemala"]),8), 
                            accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="Guatemala"]),8),
                            npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="Guatemala"]),8),
                            v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="Guatemala"]),8), 
                            v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="Guatemala"]),8),
                            rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="Guatemala"]),8),
                            impunity=rep(mean(plotdat2$impunity[plotdat2$country=="Guatemala"]),8), 
                            cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="Guatemala"]),8),
                            govteff=rep(mean(plotdat2$govteff[plotdat2$country=="Guatemala"]),8)),
                 type="response",na.action=na.pass)

#Psuedo dataset for predict loop
pm2 <- plotmodel

#Bootstrapping CI
library(mvtnorm)
library(MASS)
set.seed(5678)
mat1 <- matrix(0,nrow=8,ncol=50000) #matrix to hold predicted probability for Brazil predictions
mat2 <- matrix(0,nrow=8,ncol=50000) #matrix to hold predicted probability for Guatemala predictions
for(i in 1:50000){
  p1 <- mvrnorm(n=1,mu=plotmodel$coef,Sigma=vcov(plotmodel)) #Assume coefs are MV normally distributed with observed coefficient as the mean. 
  #Take a draw from that distribution, with mu = actual coefficient, and sigma as the vcov matrix of the model
  
  #These are now your coefficients 
  pm2$coefficients <- p1
  
  #LOOPS first for Brazil
  mat1[,i] <- predict(pm2,
                      data.frame(lamuse = lmu, 
                                 revenue_world=rep(32.36,8),
                                 postjud = factor(rep("5",8),levels=levels(factor(plotdat2$postjud))),
                                 hon=rep(0,8),
                                 per=rep(0,8),
                                 ven=rep(0,8),
                                 col=rep(0,8),
                                 gua=rep(0,8),
                                 els=rep(0,8),
                                 chl=rep(0,8),
                                 par=rep(0,8),
                                 ecu=rep(0,8),
                                 bol=rep(0,8),
                                 pan=rep(0,8),
                                 bra=rep(1,8),
                                 mex=rep(0,8),
                                 arg=rep(0,8),
                                 sur=rep(0,8),
                                 cra=rep(0,8),
                                 dmr=rep(0,8),
                                 nic=rep(0,8),
                                 tat=rep(0,8),
                                 hti=rep(0,8),
                                 pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="Brazil"]),8), 
                                 accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="Brazil"]),8),
                                 npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="Brazil"]),8),
                                 v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="Brazil"]),8), 
                                 v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="Brazil"]),8),
                                 rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="Brazil"]),8),
                                 impunity=rep(mean(plotdat2$impunity[plotdat2$country=="Brazil"]),8), 
                                 cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="Brazil"]),8),
                                 govteff=rep(mean(plotdat2$govteff[plotdat2$country=="Brazil"]),8)),
                      type="response",na.action=na.pass)
  #Now another loop, this time for Guatemala  
  mat2[,i] <- predict(pm2,
                      data.frame(lamuse = lmu, 
                                 revenue_world=rep(11.51,8),
                                 postjud = factor(rep("5",8),levels=levels(factor(plotdat2$postjud))),
                                 hon=rep(0,8),
                                 per=rep(0,8),
                                 ven=rep(0,8),
                                 col=rep(0,8),
                                 gua=rep(1,8),
                                 els=rep(0,8),
                                 chl=rep(0,8),
                                 par=rep(0,8),
                                 ecu=rep(0,8),
                                 bol=rep(0,8),
                                 pan=rep(0,8),
                                 bra=rep(0,8),
                                 mex=rep(0,8),
                                 arg=rep(0,8),
                                 sur=rep(0,8),
                                 cra=rep(0,8),
                                 dmr=rep(0,8),
                                 nic=rep(0,8),
                                 tat=rep(0,8),
                                 hti=rep(0,8),
                                 pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="Guatemala"]),8), 
                                 accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="Guatemala"]),8),
                                 npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="Guatemala"]),8),
                                 v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="Guatemala"]),8), 
                                 v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="Guatemala"]),8),
                                 rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="Guatemala"]),8),
                                 impunity=rep(mean(plotdat2$impunity[plotdat2$country=="Guatemala"]),8), 
                                 cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="Guatemala"]),8),
                                 govteff=rep(mean(plotdat2$govteff[plotdat2$country=="Guatemala"]),8)),
                      type="response",na.action=na.pass)
}
#Create bins to hold the UB and LB for the 95% CI
low1 <- rep(0,8)
low2 <- rep(0,8)
high1 <- rep(0,8)
high2 <- rep(0,8)

#Loops to get the 0.025 and 0.975 quantiles, corresponding to the bounds of the 95% CI
for(i in 1:8){
  high1[i] <- quantile(mat1[i,],0.025)
  high2[i] <- quantile(mat1[i,],0.975)
  low1[i] <- quantile(mat2[i,],0.025)
  low2[i] <- quantile(mat2[i,],0.975)
}

#PLOT - what is presented as Figure 1 in paper

dev.new(width = 330, height = 550, unit = "px")

#Panel A (Brazil)
plot(lmu,p_bra,type="l",lwd=5,xlab="Monetary award (logged)",ylab="Predicted Probability",
     ylim=c(0,0.95))
lines(lmu,high1,lty=3,lwd=4)
lines(lmu,high2,lty=3,lwd=4)

#Panel B (Guatemala)
plot(lmu,p_gua,type="l",lwd=5,xlab="Monetary award (logged)",ylab="Predicted Probability",
     ylim=c(0,0.95),col="grey50")
lines(lmu,low1,lty=3,lwd=4,col="grey50")
lines(lmu,low2,lty=3,lwd=4,col="grey50")

############################################
##TABLE 3: Results for Hypotheses 2 and 2A##
############################################

#Column 1
reg2a1 <- glm(fcomp~as.factor(postjud)+log(total_victims)+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg2a1)

#Column 2
reg2 <- glm(fcomp~as.factor(postjud)+total_victims+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
              govteff+rule_law+impunity+cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg2)

#Column 3
creg2a2 <- glm(fcomp~0+as.factor(postjud)+log(total_victims)*rep_prog+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg2a2)

#Column 4
creg2a <- glm(fcomp~0+as.factor(postjud)+log(total_victims)*rep_prog+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
                govteff+rule_law+impunity+cumj_cy_as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg2a)



######################################
##FIGURE 2: EL SALVADOR AND HONDURAS##
######################################

##log variable for model
plotdat2$ltv <- log(plotdat2$total_victims)

##same as Column 4 in Table 3
plotmodel2 <- glm(fcomp~0+as.factor(postjud)+ltv*rep_prog+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
                    govteff+rule_law+impunity+cumj_cy+hon+per+ven+col+sur+arg+cra+dmr+hti+tat+nic+
                    gua+els+chl+par+ecu+bol+pan+bra+mex,family="binomial"(link="cloglog"),data=plotdat2)

#total victims (logged)
tv <- seq(0,11,1)

#prediction model for Honduras
p_hon <- predict(plotmodel2,
                 data.frame(ltv = tv, 
                            rep_prog=rep(0,12),
                            postjud = factor(rep("5",12),levels=levels(factor(plotdat2$postjud))),
                            hon=rep(1,12),
                            per=rep(0,12),
                            ven=rep(0,12),
                            col=rep(0,12),
                            gua=rep(0,12),
                            els=rep(0,12),
                            chl=rep(0,12),
                            par=rep(0,12),
                            ecu=rep(0,12),
                            bol=rep(0,12),
                            pan=rep(0,12),
                            bra=rep(0,12),
                            mex=rep(0,12),
                            arg=rep(0,12),
                            sur=rep(0,12),
                            cra=rep(0,12),
                            dmr=rep(0,12),
                            nic=rep(0,12),
                            tat=rep(0,12),
                            hti=rep(0,12),
                            pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="Honduras"]),12), 
                            accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="Honduras"]),12),
                            npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="Honduras"]),12),
                            v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="Honduras"]),12), 
                            v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="Honduras"]),12),
                            rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="Honduras"]),12),
                            impunity=rep(mean(plotdat2$impunity[plotdat2$country=="Honduras"]),12), 
                            cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="Honduras"]),12),
                            govteff=rep(mean(plotdat2$govteff[plotdat2$country=="Honduras"]),12)),
                 type="response",na.action=na.pass)

#prediction model for El Salvador
p_els <- predict(plotmodel2,
                 data.frame(ltv = tv, 
                            rep_prog=rep(1,12),
                            postjud = factor(rep("5",12),levels=levels(factor(plotdat2$postjud))),
                            hon=rep(0,12),
                            per=rep(0,12),
                            ven=rep(0,12),
                            col=rep(0,12),
                            gua=rep(0,12),
                            els=rep(1,12),
                            chl=rep(0,12),
                            par=rep(0,12),
                            ecu=rep(0,12),
                            bol=rep(0,12),
                            pan=rep(0,12),
                            bra=rep(0,12),
                            mex=rep(0,12),
                            arg=rep(0,12),
                            sur=rep(0,12),
                            cra=rep(0,12),
                            dmr=rep(0,12),
                            nic=rep(0,12),
                            tat=rep(0,12),
                            hti=rep(0,12),
                            pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="El Salvador"]),12), 
                            accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="El Salvador"]),12),
                            npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="El Salvador"]),12),
                            v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="El Salvador"]),12), 
                            v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="El Salvador"]),12),
                            rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="El Salvador"]),12),
                            impunity=rep(mean(plotdat2$impunity[plotdat2$country=="El Salvador"]),12), 
                            cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="El Salvador"]),12),
                            govteff=rep(mean(plotdat2$govteff[plotdat2$country=="El Salvador"]),12)),
                 type="response",na.action=na.pass)

#Psuedo dataset for predict loop
pm3 <- plotmodel2

#Bootstrapping CI
set.seed(5678)
mat3 <- matrix(0,nrow=12,ncol=50000) #matrix to hold predicted probability for El Salvador
mat4 <- matrix(0,nrow=12,ncol=50000) #matrix to hold predicted probability for Honduras
for(i in 1:50000){
  p1 <- mvrnorm(n=1,mu=plotmodel2$coef,Sigma=vcov(plotmodel2)) #Assume coefs are MV normally distributed with observed coefficient as the mean. 
  #Take a draw from that distribution, with mu = actual coefficient, and sigma as the vcov matrix of the model
  
  #These are now your coefficients 
  pm3$coefficients <- p1
  
  #LOOPS - predict first for El Salvador
  mat3[,i] <- predict(pm3,
                      data.frame(ltv = tv, 
                                 rep_prog=rep(1,12),
                                 postjud = factor(rep("5",12),levels=levels(factor(plotdat2$postjud))),
                                 hon=rep(0,12),
                                 per=rep(0,12),
                                 ven=rep(0,12),
                                 col=rep(0,12),
                                 gua=rep(0,12),
                                 els=rep(1,12),
                                 chl=rep(0,12),
                                 par=rep(0,12),
                                 ecu=rep(0,12),
                                 bol=rep(0,12),
                                 pan=rep(0,12),
                                 bra=rep(0,12),
                                 mex=rep(0,12),
                                 arg=rep(0,12),
                                 sur=rep(0,12),
                                 cra=rep(0,12),
                                 dmr=rep(0,12),
                                 nic=rep(0,12),
                                 tat=rep(0,12),
                                 hti=rep(0,12),
                                 pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="El Salvador"]),12), 
                                 accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="El Salvador"]),12),
                                 npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="El Salvador"]),12),
                                 v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="El Salvador"]),12), 
                                 v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="El Salvador"]),12),
                                 rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="El Salvador"]),12),
                                 impunity=rep(mean(plotdat2$impunity[plotdat2$country=="El Salvador"]),12), 
                                 cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="El Salvador"]),12),
                                 govteff=rep(mean(plotdat2$govteff[plotdat2$country=="El Salvador"]),12)),
                      type="response",na.action=na.pass)
  #prediction for Honduras
  mat4[,i] <- predict(pm3,
                      data.frame(ltv = tv, 
                                 rep_prog=rep(0,12),
                                 postjud = factor(rep("5",12),levels=levels(factor(plotdat2$postjud))),
                                 hon=rep(1,12),
                                 per=rep(0,12),
                                 ven=rep(0,12),
                                 col=rep(0,12),
                                 gua=rep(0,12),
                                 els=rep(0,12),
                                 chl=rep(0,12),
                                 par=rep(0,12),
                                 ecu=rep(0,12),
                                 bol=rep(0,12),
                                 pan=rep(0,12),
                                 bra=rep(0,12),
                                 mex=rep(0,12),
                                 arg=rep(0,12),
                                 sur=rep(0,12),
                                 cra=rep(0,12),
                                 dmr=rep(0,12),
                                 nic=rep(0,12),
                                 tat=rep(0,12),
                                 hti=rep(0,12),
                                 pobj_count=rep(mean(plotdat2$pobj_count[plotdat2$country=="Honduras"]),12), 
                                 accept_resp=rep(mean(plotdat2$accept_resp[plotdat2$country=="Honduras"]),12),
                                 npecuniary=rep(mean(plotdat2$npecuniary[plotdat2$country=="Honduras"]),12),
                                 v_achr4=rep(mean(plotdat2$v_achr4[plotdat2$country=="Honduras"]),12), 
                                 v_achr5=rep(mean(plotdat2$v_achr5[plotdat2$country=="Honduras"]),12),
                                 rule_law=rep(mean(plotdat2$rule_law[plotdat2$country=="Honduras"]),12),
                                 impunity=rep(mean(plotdat2$impunity[plotdat2$country=="Honduras"]),12), 
                                 cumj_cy=rep(mean(plotdat2$cumj_cy[plotdat2$country=="Honduras"]),12),
                                 govteff=rep(mean(plotdat2$govteff[plotdat2$country=="Honduras"]),12)),
                      type="response",na.action=na.pass)
}
#Create bins to hold the UB and LB for the 95% CI
low3 <- rep(0,12)
low4 <- rep(0,12)
high3 <- rep(0,12)
high4 <- rep(0,12)

#Loops to get the 0.025 and 0.975 quantiles, corresponding to the bounds of the 95% CI
for(i in 1:12){
  high3[i] <- quantile(mat3[i,],0.025)
  high4[i] <- quantile(mat3[i,],0.975)
  low3[i] <- quantile(mat4[i,],0.025)
  low4[i] <- quantile(mat4[i,],0.975)
}

#PLOT - Figure 2

dev.new(width = 330, height = 550, unit = "px")

##Panel A
plot(tv,p_els,type="l",lwd=5,xlab="Total victims (logged)",ylab="Predicted Probability",
     ylim=c(0,0.3))
lines(tv,high3,lty=3,lwd=4)
lines(tv,high4,lty=3,lwd=4)

##Panel B
plot(tv,p_hon,type="l",lwd=5,xlab="Total victims (logged)",ylab="Predicted Probability",
     ylim=c(0,0.3),col="grey50")
lines(tv,low3,lty=3,lwd=4,col="grey50")
lines(tv,low4,lty=3,lwd=4,col="grey50")

#########################
##TABLE 4: Hypothesis 3##
#########################

#Column 1
reg3a1 <- glm(fcomp~0+as.factor(postjud)+type97+type98+type99+as.factor(country),
              data=cd5[cd5$pnp==1 & cd5$h3test==1,],family="binomial"(link="cloglog"))
summary(reg3a1)

#Column 2
reg3 <- glm(fcomp~0+as.factor(postjud)+type97+type98+type99+log(amount_use)+log(total_victims)+
              pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
              govteff+rule_law+impunity+cumj_cy+as.factor(country),
            data=cd5[cd5$pnp==1 & cd5$h3test==1,],family="binomial"(link="cloglog"))
summary(reg3)

#Column 3
creg3a <- glm(fcomp~0+as.factor(postjud)+type97+type98+type99+log(amount_use)+log(total_victims)+
                rep_prog+revenue_world+pobj_count+accept_resp+
                npecuniary+v_achr4+v_achr5+govteff+rule_law+impunity+cumj_cy+as.factor(country),
              data=cd5[cd5$pnp==1 & cd5$h3test==1,],family="binomial"(link="cloglog"))
summary(creg3a)

#Column 4
creg3ba2 <- glm(fcomp~0+as.factor(postjud)+type97+type98+type99+type2a+type3a+
                  as.factor(country),
                data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg3ba2)

#Column 5
creg3b <- glm(fcomp~0+as.factor(postjud)+type97+type98+type99+type2a+type3a+
                log(amount_use)+log(total_victims)+rep_prog+revenue_world+pobj_count+accept_resp+
                npecuniary+v_achr4+v_achr5+govteff+rule_law+impunity+cumj_cy+as.factor(country),
              data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg3b)


######################################
##APPENDIX TABLE A.1 (Omnibus Model)##
######################################

#Column 1
reg4a <- glm(fcomp~0+as.factor(postjud)+log(amount_use)+log(total_victims)+
               pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+govteff+rule_law+impunity+
               cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg4a)

#Column 2
reg4b <- glm(fcomp~0+as.factor(postjud)+log(amount_use)+log(total_victims)+revenue_world+rep_prog+
               pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+govteff+rule_law+impunity+
               cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg4b)

#Column 3
reg4c <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*revenue_world+log(total_victims)+rep_prog+
               pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
               govteff+rule_law+impunity+cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg4c)

#Column 4
reg4d <- glm(fcomp~0+as.factor(postjud)+log(total_victims)*rep_prog+log(amount_use)+revenue_world+
               pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
               govteff+rule_law+impunity+cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(reg4d)



######################################################
##APPENDIX TABLE A.2 (Alternative capacity measures)##
######################################################

#Column 1
creg1b1 <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*gdp_cap_cur+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg1b1)

#Column 2
creg1b <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*gdp_cap_cur+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
                impunity+cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg1b)

#Column 3
creg1c1 <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*gdppc_q2a+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg1c1)

#Column 4
creg1c <- glm(fcomp~0+as.factor(postjud)+log(amount_use)*gdppc_q2a+pobj_count+accept_resp+npecuniary+v_achr4+v_achr5+
                impunity+cumj_cy+as.factor(country),data=cd5[cd5$pnp==1,],family="binomial"(link="cloglog"))
summary(creg1c)


